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We model the dynamics of avalanches in granular assemblies in partly filled rotating cylinders 
using a mean-field approach. We show that, upon varying the cylinder angular velocity uj, the 
system undergoes a hysteresis cycle between an intermittent and a continuous fiow regimes. In the 
intermittent fiow regime, and approaching the transition, the avalanche duration exhibits critical 
slowing down with a temporal power-law divergence. Upon adding a white noise term, and close to 
the transition, the distribution of avalanche durations is also a power-law. The hysteresis, as well 
as the statistics of avalanche durations, are in good qualitative agreement with recent experiments 
in partly filled rotating cylinders. 
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Introduction. — The dynamics of granular materials 
is not well understood in spite of its widespread scien- 
tific and technological interest Bistability, segrega- 
tion, arching, hysteresis, instabilities and other proper- 
ties found in noncohesive granular assemblies make these 
systems very difBcult to analyze. For over two centuries, 
most studies have focused on the two extreme regimes: 
the fluidlike continuous flow and the static compact state. 
The intermediate quasistatic regime, with intermittent 
transitions between flowing and static behavior has re- 
cently generated a flurry of activity including several neat 
experiments [ij- [|| and novel theoretical proposals (see, 
e.g., Ref. |§). 

Several experiments ^ have studied the dynam- 
ics of granular assemblies in partly filled rotating drums. 
In particular, Rajchenbach monitored the following 
quantities: (z) the statistics of avalanches as a function of 
the rotating speed of the cylinder, uj, (ii) the transition 
from the discrete to the steady regime, including exper- 
imental evidence of hysteresis as a function of w , and 
(in) the dependence of the current of particles, j, with 
the angle 9 between the top layer and the horizontal. It 
is the purpose of this work to study these issues by using 
a simple mean-field model with two dynamical variables. 
The recent studies on this subject raise many interesting 
questions in granular transport, like hysteresis in w (see, 
for instance, fig. 1 of ref. Q), and avalanche statistics, 
and this work addresses some of them. 

Our mean-field approach reduces the many degrees 
of freedom of a granular material in a rotating drum 
into two: the average velocity, v (= dx/dt = i), of 
the avalanche and the average angle, 0, of the granu- 
lar assembly top surface (of length L, < x < L). 
The grains move in the "downhill" x-direction of motion, 
which forms an angle 9 with the horizontal, and the shear 
between two layers with a separation d is 7 = v/d. An 
important feature of our approach is that the dynamics of 
the model used here predicts the occurrence of avalanches 
as a function of lu. By slowly increasing to, we obtain 
the angular velocity, ijJ+, at which the system undergoes 
a transition from the discrete flow regime to a contin- 
uous flow (see fig. 1). At this point, if lu is decreased, 
the reverse transition is observed at , where U- < uj+ . 
Thus, our model exhibits hysteresis in lo during the tran- 
sition between the "intermittent" and "continuous" flow 
regimes. 

Bagnold extensively studied the dissipation pro- 
duced by the grain collisions in several regimes, including 
the quasistatic intermittent state, finding in the latter a 
periodic response over time. On the other hand, sim- 
ple rules for the dynamics |^ do not predict a periodic 
response over time but the absence of length and time 
scales. Our model can produce distributions in agree- 
ment with recent experiments [^Q. Furthermore, when 
u —f uj^ , we predict that the system exhibits a power-law 
critical slowing down. 



Dissipation in Granular Media. — In granular fiow, 
the friction force depends on v^ for large values of the 
velocity. The models of Bagnold ||] and Jenkins and 
Savage ||l^ correctly describe this limit. Jaeger et al 
analyzed how the stress varies as v approaches zero, 
obtaining the following form for the total friction force: 
mv^/2A -I- amgcos9/K where k = 1 + aiv'^ / gdcos9 ; 
which in dimensionless form can be written as: F = 
fiv^ + acos9/ K, , where now k = 1 + aiv^ / cos9, v — 



v/y/dg , and a, ai, (3, X are constants. Here, the first 
(second) term dominates the friction force at high (low) 
velocities, and, when u — s- 0, the second term represents 
the maximum dry friction due to rubbing forces. 

Dynamics. — The equation of motion is v = a — F, 
where F = F{9, v) is the friction force, and a is the grav- 
itational force ((X sin 6*). If c = aai/P > 1, then 
the frictional force decreases as a function of u, until a 
minimum is reached at Vr, and then increases. We fo- 
cus on this case, c > 1, illustrated in fig. 2(a). The 
stationary points of F determine the maximum, 9m, and 
minimum (or repose) 9r, angle of stability. Explicitly: 
9m = arctan(Q;), and 9r = arctan(a(2Y^ — l)/c). 

Let us consider a stationary {v = 0) granular assem- 
bly with a horizontal {9 = 0) top surface contained in 
a partly filled cylinder. If we slowly rotate the cylinder 
at an angular velocity 9 = lu, the average angle of the 
granular assembly surface, 9 grows until it reaches the 
maximum angle of stability, 9m, where the system sud- 
denly produces avalanches. During this sudden transport 
of mass, the average velocity of the sand grains, v, be- 
comes nonzero and 9 decreases until it reaches the angle 
of repose, 9r, where the system relaxes {v=0). At this 
point the cycle repeats itself. 

In order to model this dynamics, it is useful to know 
the time evolution of the mean angle 9. For this pur- 
pose, let us consider the two extreme heights (shown in 
fig. 2(b)) of granular vertical columns. Since the column 
heights are proportional to the mass of the columns, the 
rate of mass transport is proportional to the hcigth dif- 
ference: {0 < X < L), 
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The proportionality constant is the inverse of the prop- 
agation time. The sum of these equations give the mass 
conservation condition dt{hi + ft.2) = 0; which is valid in 
a closed rotating drum, while substracting them gives 



= -2— tan( 

1j 



(2) 



since /ii — /12 = ^ sin 9. 

In this work, we concentrate our efforts in this dx9 ^ 
case. Elsewhere, we will consider the case with local 
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spatial variations in 9{x,t), relevant to describe the S- 
shaped curves seen in the continuous flow regime for very 
large values of u. 

We thus model the dynamics of avalanches with the 
following equations of motion: 

i) ^ sine-l3v^ '^^^^ = s\ne-F(e,v) (3) 

1 + aiv-^ / cosy 

9 = -pv tan 6 + + ri{t) (4) 

where p = eq.(3) is Newton's equation of mo- 

tion, and the very small Gaussian noise term -qit) simu- 
lates small fluctuations in cu and satisfles: < ri{t) >— 
and < r]{t)r]{t') Dq S{t -t') , with Dq < uj. 

In the continuous flow regime, the equation for the 
steady state (w = 0, = 0) is 

sin 6* = /3 I — I -I- a cos^ 6* ( cos 6* -I- ai ( — ) | 

(5) 

This solution corresponds to any point in the branch DC' 
in fig. 2(a) with positive slope in F{v). Let us denote by 
Vr (vm) the velocity corresponding to 9r {9m) in eq. (3) 
when i) = 0. When lo <C pvr tan 9r , the system has 
intermittent avalanches, and in each one of them eq.(4) 
predicts that 9 relaxes fast to 9r- When lo > pVm tan0„i , 
the system is certainly in the continuous flow regime at 
some point C" above C (see fig. 2(a)). If pvrta,n9r < 
UJ < pw„jtan0„j, then the system evolves from C to a 
point between C and D and either remains there (if the 
stationary solution is stable at this point) or jumps back 
(i.e., 9 ^ 9r and u ^ 0) to the AB segment and the 
avalanche stops. If w < purtan^r, then the system 
certainly decays to the AB segment and has intermittent 
avalanches evolving through the loop ABCD, and so on, 
since D is unstable. 

Numerical Results. — Figure 1 shows the time depen- 
dence of: (a) the flow velocity v, (b) the angle 9 of the 
top layer, and (c) the angular velocity to of the rotat- 
ing drum, obtained by numerically solving eqs.(3-4) with 
Dq = 0. For each avalanche, v has a peak and 9 has a 
sudden decrease from 9m ~ 0.8 to 9r ~ 0.65. The v = 
regions in between peaks correspond to the loading time 
in between avalanches, where the angle increases linearly 
in time from 9r to 9m- For clarity, only three avalanches 
are shown. They are approaching the transition point, 
cij+, to the continuous flow regime, and their periods be- 
come longer close to lu^ , i. e. the system exhibits a critical 
slowing down. The values of both a;+ and lu^ depend on 
Lu, and the constants a, ai, /?, and Dq in eqs.(3-4). The 
addition of noise also reduces the hysteresis. 

Notice the asymmetric form of the velocity peaks in 
flg. 1(a). During a loading time, ta, the system moves 
from A to B, on the branch AB in fig. 2(a). Right after- 
wards, an avalanche starts and v has: a sudden increase 



lasting a time tb {B ^ C in fig. 2 (a)), a slow decrease 
with duration tc {C ^ D), and a final decrease lasting 
a time tjj {D ^ A). 

Equations (3-4), without a random source term, and 
for a fixed w, produce a periodic sequence of avalanches. 
However, the real many-body system inevitably has dis- 
order. When Dq ^ 0, we observe a richer behavior (see 
figs. 3 and 4(b)) with broad distributions. 

The statistics of the avalanche durations are shown in 
Fig. 3. In order to gain more insight into the several time 
scales involved in the problem, we show the distributions, 
N, for (a) TA, (b) tb, (c) tc, and (d) td, for 3x10^'' 
and Dq ^ 10"*. 

Figure 4(a) shows the evolution of td when lu — s- 
UJ+. We have considered two different possibilities for 
the divergence: logarithmic (left vertical axis of the inset) 
and power-law (right vertical axis). For Dq = 0, the 
behavior of tb at the transition is found to be 

Tb ^ (w+ — ^)^^ for ^ < ^+ ■ (6) 

Temporal critical slowing down behavior has also been 
found in bistable optical systems ||ll|]. From eq.(6) and 
recalling that + 77 is a Gaussian random variable, one 
obtains that N{tb) ~ t^^. This is checked numerically 
in fig. 4(b) where we have obtained the distribution of tb 
for fixed uj = 4.6 x 10"'^, which is slightly below uj^{Dq = 
0), with Dq = 10"*. 

In the intermittent flow regime, our dynamical model 
exhibits an attractive limit cycle around an unstable fixed 
point. For 9 < 9m , v = is a, stable attractor manifold. 
In the continuum flow regime, the dynamical system is 
in an attractive fixed point with nonzero velocity. Here, 
a harmonic approximation to Newton's equation (eq.(3)) 
with V — gives 

v-Vr-{9- 9r)^l^ (7) 

for the relationship between the mass current flow and 
the angle, for oj > uj-. This result is consistent with 
experimental results by Rajchenbach [||. 

Limitations. — Our mean-field equations have only two 
degrees of freedom for modeling an extremely complex 
system. Thus, it might need to be extended in order to 
provide a more complete description of the system. The 
beauty of the model lies precisely in its simplicity and in 
the fact that it can naturally describe several dynamical 
features observed in experiments in granular media. 

Conclusions. — We have studied a mean- field model of 
granular assemblies driven to the threshold of instability, 
where they produce avalanches. From this simple model, 
we obtained the distributions of durations of avalanches, 
the transition between the two flowing regimes, and the 
hysteresis between them in a natural manner, and with a 
good qualitative agreement with results recently obtained 
using experiments in rotating cylinders. Furthermore, 
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when Lo lj^, the avalanche durations exhibit a critical 
slowing down with a temporal power-law divergence. 
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FIG. 2. (a) Schematic plot of the friction force versus ve- 
locity V for the case c > 1. An initially flat (6 = 0) and 
stationary {v = 0) granular surface is tilted with an angular 
velocity uj, thus having the following sequence of values for 
the friction force: O ^ A —f B, from the F — v — origin, 
O. The velocity, v, is zero in the stick-slip case with infinite 
slope for the OAB segment. In the AB branch, the loading 
time, TA, satisfies < ta >= {Om — Qr)l^. When the maximum 
angle of stability is reached in B, an avalanche {v ^ 0) starts 
and the system jumps to C. Different situations, described 
in the text, can now develop according to the value of oj. (b) 
Schematic diagram of a rotating drum indicating the extreme 
heigths, the mean downhill velocity v, and the mean angle B. 

FIG. 3. Statistics of the avalanche durations. Distri- 
butions, iV, for (a) TA, (b) tb, (c) tc, and (d) to, for 
a; = 3 X IQ-* and Do = 10"^ 

FIG. 4. (a) Time evolution of td(cij) 
when LO — > (i)+(_Do = 0). This critical slowing down with 
a power law divergence is shown in the inset by the dashed 
line corresponding to the right vertical logarithmic axis. For 
comparison purposes, the solid line in the inset, corresponding 
to the left vertical linear axis, has a small curvature because 
of the small value of the power (p ~ 0.34). (b) Distribution 
of To for = 4.6 X 10"'', which is slightly below 01+ (Do = 0). 
The inset shows the same data in semi-log and log-log scales. 



FIG. 1. Shows the time dependence of: (a) the flow velocity 
V, (b) the angle Q of the top layer, and (c) the angular velocity 

of the rotating drum. Throughout this work, we use the 
parameters values: a = 1, ai = 2, /? = 1/2, and p = 10~^. 
For these ramping rates, we obtain a;+ = 4.76 x 10"'' and 
LO- = 4.66 X lO^'', for the Do = deterministic case. 
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